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ABSTRACT 

We  develop  and  evaluate  a  suite  of  advanced  algorithms  which  provide  significantly-improved  capabilities  for  finding, 
fixing,  and  tracking  multiple  ballistic  and  flying  low  observable  objects  in  highly  stressing  cluttered  environments. 
The  algorithms  have  been  developed  for  use  in  satellite-based  staring  and  scanning  optical  surveillance  suites  for  ap¬ 
plications  including  theatre  and  intercontinental  ballistic  missile  early  warning,  trajectory  prediction,  and  multi-sensor 
track  handoff  for  midcourse  discrimination  and  intercept.  The  system  performs  sensor  motion  compensation  providing 
sub-pixel  stabilization  (to  1/100  of  a  pixel),  as  well  as  advanced  spatial-temporal  clutter  estimation  and  suppression  to 
below  sensor  noise  levels,  followed  by  statistical  background  modeling  and  nonlinear  Bayesian  multiple-target  track- 
before-detect  filtering.  The  multiple-target  tracking  is  performed  in  physical  world  coordinates  to  allow  for  multi¬ 
sensor  fusion,  trajectory  prediction,  and  intercept.  The  algorithms  are  designed  to  handle  a  wide  variety  of  real-world 
challenges:  the  scene  background  may  contain  significant  celestial,  earth  limb,  or  terrestrial  clutter,  and  the  targets  of 
interest  may  be  dim,  relative  to  the  scene  background.  The  performance  of  the  developed  algorithms  is  demonstrated 
using  real-world  data  containing  resident  space  objects  observed  from  the  MSX  platform,  with  backgrounds  varying 
from  celestial  to  combined  celestial  and  earth  limb.  Simulation  results  are  also  presented  for  parameterized  varia¬ 
tions  in  signal-to-clutter  levels  (down  to  1/1000)  and  signal-to-noise  levels  (down  to  1/6)  for  simulated  targets  against 
real-world  terrestrial  clutter  backgrounds. 


1.  INTRODUCTION 

The  problem  of  efficient  detection  and  tracking  of  multiple  dim  targets  is  a  challenge  for  space-based  IR/EO  sensors. 
In  fact,  in  the  raw  image  sequence  the  signal-to-clutter  ratio  (SCR)  may  be  as  low  as  10-3,  in  which  case  detection 
and  tracking  without  preliminary  spatial-temporal  image  processing  and  clutter  rejection  is  almost  impossible.  For 
geostationary  sensors,  this  problem  has  been  partially  solved  by  Tartakovsky  and  Brown  [4].  As  has  been  shown 
in  this  work,  the  excess  false  alarm  problem  can  be  addressed  only  through  sophisticated  spatial-temporal  image 
processing  for  clutter  rejection  (CLUR).  Specifically,  conventional  spatial  and  differencing  methods  are  not  sufficient, 
and  advanced  spatial-temporal  CLUR  filters  are  required.  One  of  the  major  nuisance  factors  (even  for  mechanically 
stabilized  IR/EO  sensors  and  geostationary  platforms)  that  prevents  efficient  temporal  processing  is  LOS  (line-of- 
sight)  non-regular  and  uncontrollable  changes  due  to  vibration,  which  results  in  translational  and  rotational  image 
distortions.  Existing  mechanical  methods  of  sensor  stabilization  are  very  expensive  and  are  not  efficient  enough,  as 
was  demonstrated  in  [4] . 

This  paper  focuses  on  the  development  and  evaluation  of  a  suite  of  advanced  algorithms  which  provide  significantly- 
improved  capabilities  for  finding,  fixing,  and  tracking  multiple  ballistic  and  flying  low  observable  objects  in  highly 
stressing  cluttered  environments.  The  algorithms  have  been  developed  for  use  in  satellite-based  staring  and  scanning 
optical  surveillance  suites  for  applications  including  theatre  and  intercontinental  ballistic  missile  early  warning,  trajec¬ 
tory  prediction,  and  multi-sensor  track  handoff  for  midcourse  discrimination  and  intercept.  As  illustrated  in  Figure  1 , 
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the  functions  performed  by  the  algorithms  include  sensor  motion  (and/or  platform  vibration)  compensation,  providing 
sub-pixel  stabilization  (to  1/100  of  a  pixel),  advanced  spatial-temporal  clutter  estimation  and  suppression  to  below  sen¬ 
sor  noise  levels,  statistical  background  modeling,  and  Bayesian  nonlinear  multiple-target  track-before-detect  filtering. 
The  multiple-target  tracking  is  performed  in  physical  world  coordinates  to  allow  for  multi-sensor  fusion,  trajectory 
prediction,  and  intercept.  Output  of  detected  object  cues  and  data  visualization  are  also  provided. 


Figure  1.  Block  diagram  illustrating  functionality  of  the  enhanced  low  observables  algorithms. 

The  algorithms  are  designed  to  handle  a  wide  variety  of  real-world  challenges.  Sequences  of  images  may  be  highly 
complex  and  infinitely  varied — the  scene  background  may  contain  significant  celestial,  earth  limb,  or  terrestrial  clutter. 
For  example,  when  viewing  combined  earth  limb  and  terrestrial  scenes,  a  combination  of  stationary  and  non- stationary 
clutter  may  be  present,  including  cloud  formations,  varying  atmospheric  transmittance  and  reflectance  of  sunlight  and 
other  celestial  light  sources,  aurora,  glint  off  sea  surfaces,  and  varied  natural  and  man-made  terrain  features.  The 
targets  of  interest  may  be  (and  usually  are)  dim,  relative  to  the  scene  background,  rendering  much  of  the  existing 
deployed  software  useless  for  optical  target  detection  and  tracking.  Additionally,  it  may  be  necessary  to  detect  and 
track  a  large  number  of  objects  in  the  threat  cloud,  and  these  objects  may  not  always  be  resolvable  in  individual  data 
frames. 

The  performance  of  the  developed  algorithms  is  demonstrated  using  real-world  data  containing  resident  space 
objects  observed  from  the  MSX  platform,  with  backgrounds  varying  from  celestial  to  combined  celestial  and  earth 
limb.  Simulation  results  are  also  presented  for  parameterized  variations  in  signal-to-clutter  levels  (down  to  1/1000) 
and  signal-to-noise  levels  (down  to  1/6)  for  simulated  targets  against  real-world  terrestrial  clutter  backgrounds.  We 
also  discuss  algorithm  processing  requirements  and  software  processing  capabilities  from  our  on-going  development 
of  an  image  processing  toolkit  (iPTK). 

The  rest  of  the  paper  is  organized  as  follows.  In  Section  2,  we  outline  the  developed  CLUR  methods  for  geostation¬ 
ary  platforms,  and  illustrate  that  this  development  is  indeed  important  by  comparing  with  the  typical  state-of-practice 
differencing  methods.  Section  3  is  devoted  to  the  clutter  suppression  issues  for  highly  non- stationary  conditions  that 
arise  in  scenarios  where  sensors  are  situated  on  moving  platforms  such  as  low  Earth  orbit  (LEO)  satellites,  highly 
elliptical  orbit  (HEO)  satellites  (high  Earth  orbits),  etc.  In  Section  4,  we  outline  our  non-linear  track-before-detect 
algorithm,  and  we  provide  example  results  in  stressing  cases,  including  tracking  objects  with  dramatically  different 
velocities  and  intensities  in  MSX  SB  V  data,  and  3D  coordinate  tracking  of  an  extremely  dim  simulated  ballistic  target 
against  a  EUMETSAT  clutter  background.  Finally,  in  Section  5  we  draw  conclusions  and  make  final  remarks. 

2.  CLUTTER  REJECTION  FOR  GEOSTATIONARY  SENSORS 

In  a  geostationary  scenario,  the  developed  baseline  CLUR  technique  is  based  on  a  multi-parametric  approximation  of 
clutter  that  leads  to  an  adaptive  spatial-temporal  filter.  The  results  of  an  extensive  experimental  study  [4]  show  that  the 
proposed  algorithm  gives  a  substantial  gain  compared  to  the  best  existing  spatial  techniques  as  well  as  to  the  industry 
standard  temporal  differencing  method. 


We  observe  a  sequence  of  Nx  x  Ny  IR  images  (intensities)  that  have  the  form 

Kn 

Ai(fc)S(rij  -  rn(fe)  -  $n)  +  6„(rii  -  <S„)  +^„(rij),  n  =  1,2, . .. ,  (1) 

fe=i 

where  £n(r^-)  is  sensor  noise;  fr(r^  )  is  clutter  (background);  An(k)S(rij  —  r n(k))  is  a  signal  from  the  k- th  target  with 
spatial  coordinates  rn(fc)  =  (Xn(k),Yn(k))  and  maximal  intensity  An(k)\  S'(r^)  is  the  target  signature  related  to  the 
sensor’s  point  spread  function  (PSF);  Kn  is  an  unknown  number  of  targets  in  the  n-th  frame;  Sn  =  (5x(ri),5y(ri))  is 
an  unknown  2D  shift  due  to  the  jitter  (LOS  vibrations);  r^-  =  (xi,yj)  is  the  pixel  in  the  plain  image  with  coordinates 

{%i 1 5  •  •  •  5  Nx ,  j  1 , . . . ,  Ny . 

We  will  use  boldface  for  a  vector/matrix  notation,  e.g.,  In  =  {/n(r^),  i  =  1, . . . ,  j  =  1, . . . ,  Ny}. 

The  goal  is  to  build  a  spatial-temporal  filter  that  rejects  clutter  as  much  as  possible  and  simultaneously  compensates 
for  the  jitter  (stabilizes  the  scene),  while  preserving  signals  from  targets  as  much  as  possible. 

We  consider  residual-type  filtering  algorithms  of  the  form  /n(r^)  =  /n(r^-)  —  6n(r^ ),  where  /n(r^-)  is  the  output 
of  the  filter  and  bn(^ij)  is  an  estimate  of  clutter  6n(r^)  in  the  pixel  (i,  j)  of  the  current  n-th  frame  based  on  the 
previous  frames  Jn_m+i(r^  ), . . . ,  Jn(r^-)  in  a  sliding  window  of  the  size  m  (n  >  m). 

We  assume  that  in  the  time  interval  of  m  frames  the  function  bn( r)  is  slowly  varying,  i.e.,  we  can  neglect  changes 
due  to  physical  causes  such  as  wind,  illumination  conditions,  temperature,  convection,  etc.  Under  this  assumption,  the 
clutter  function  6n(r^)  changes  in  time  in  each  pixel  only  because  of  the  sensor  vibrations  (locally  within  a  certain 
time  interval). 

The  proposed  CLUR  algorithms  are  adaptive  and  robust  (use  minimum  prior  information  regarding  the  statisti¬ 
cal  properties  of  clutter  and  jitter).  We  suppose  that  translations  Sn  are  arbitrary  unknown  variables  bounded  by  a 
maximum  possible  amplitude  of  the  sensor  vibrations  and  6n(r)  is  an  arbitrary  unknown  non-negative  function  with  a 
bounded  spatial  band,  slowly  changing  in  time.  In  the  geostationary  staring  IR  sensor  scenario  the  value  of  Sn(r)  =  Sn 
is  just  the  parallel  shift  that  does  not  depend  on  r. 

The  value  of  m  defines  the  length  of  the  interval  in  which  clutter  does  not  change  substantially  (slowly  varying). 
Spatial-temporal  filtering  for  clutter  rejection  is  based  on  the  data  (In_m+i, . . . ,  In)  in  the  sliding  window  m. 

2.1.  Clutter  Suppression.  We  now  proceed  with  a  description  of  the  basic  idea  and  a  generic  CLUR  algorithm. 
Assume  that  the  function  bn  (r)  can  be  approximated  by  a  parametric  model 

M 

Mr)  ~  (2) 

k= i 

where  Ok  are  unknown  parameters  and  fk  (r)  are  given  functions  chosen  from  the  best  fitting  criterion. 

Let  6k(n)  denote  an  estimate  of  Ok  based  on  the  data  In_m+i, . . . ,  In  in  the  time  window  \n  —  m  +  1,  n\  of  the 
length  m.  According  to  (2),  for  any  shift  d ,  the  prediction  estimate  of  the  background  at  time  n  has  the  form: 

M 

bn(r  -  6)  =  ^2§k(n)fk(r  -  S).  (3) 

k=  1 


Assume  that  the  shifts  Ss  are  fixed  (estimated)  and  let  Ss  denote  these  estimates.  Write 

Nx  Ny  n  M  2 

£n(d,  {5S})  =  E  E  E  (Ar«)  -  E  -  U)  • 

i=  1  j= 1  s=n— ra+1  k= 1 

The  estimate  0n  is  found  from  the  least  squares  criterion: 

0n  =  argmin£n(0,  {£s}).  (4) 

e 

The  parameter  estimation  algorithm  requires  a  reasonably  accurate  estimation  of  the  shift  d ,  i.e.,  jitter  compensa¬ 
tion.  This  latter  estimation/compensation  can  be  done  either  by  an  independent  jitter  estimation  algorithm  or  iteratively 
in  the  course  of  estimating  parameters  0. 


Suppose  now  that  the  estimate  bn_i  (r)  of  the  (n  —  l)-st  frame  is  already  obtained.  For  jitter  estimation  in  the  n-th 
frame,  we  use  the  minimum  distance  estimate  fin  =  ( Sx(n ),  Sy(n )),  which  is  the  solution  of  the  following  nonlinear 
minimization  problem: 


Nx  Ny 


argmin 

|<5|^<5max  j=X 


(5) 


We  now  summarize  the  generic  CLUR  method  that  uses  an  iterative  procedure  for  jitter  compensation  (frame 
alignment)  and  the  background  estimation.  The  algorithm  starts  with  initialization,  which  provides  the  pilot  estimates 
6k(m ),  (<5i, . . . ,  fim),  and  6m(r^)  =  Y^k= 1  —  Sm).  Initialization  schemes  include  autonomous  algo¬ 

rithms  of  estimation  of  shifts  between  two  frames  based  on  simple  spline-approximations.  After  initialization  a  typical 
step  of  the  algorithm  includes  the  following  operations: 

1.  Jitter  Estimation.  The  estimate  frn-i(Uj)  obtained  from  the  previous  step  is  compared  with  the  n-th  frame  and 
the  estimate  of  jitter  8n  is  computed  as  the  solution  of  the  nonlinear  optimization  problem  (5)  with  6n_i(r^  —  (5)  = 

Efcii  Ok(n  -  l)/fc(r<j  -  5). 

2.  Estimation  of  Parameters .  Having  the  estimates  5n_m+i, . . . ,  8n,  the  estimates  6k{n)  are  computed  for  the 
n-th  frame  from  the  minimization  problem  (4).  This  recomputing  in  the  corrected  coordinate  system  is  equivalent  to 
frame  alignment. 

3.  Clutter  Estimation.  Using  the  estimates  obtained  from  (4)  and  (5),  compute  the  estimate  of  6n(r^)  for  all  i  and 
j  in  the  corrected  coordinate  system,  i.e.,  6n(r^)  =  Y^k= i  @k{n)fk{rij  ~  fin)- 

4.  Clutter  Rejection.  Using  the  estimate  6n(r^)  from  the  previous  step,  compute  the  residuals  (filtered  back¬ 
ground) 

M 

■)  fin)  =  Inijij)  ^  ^  @k  (jt)fk  (rij  fin)- 

k= 1 

Different  filters  in  the  bank  of  CLUR  filters  correspond  to  different  choices  of  the  basis  /.  For  example,  we  used 
Fourier  basis,  wavelets,  and  splines  (cubic  and  bi-linear),  as  well  as  certain  polynomial  approximations.  For  further 
details  we  refer  to  [4]. 


2.2.  Detection  and  Tracking.  The  processed  data  from  the  CLUR  filter  come  to  the  Target  Detection  Block.  The 
output  of  this  block  are  the  instantaneous  detections  (blips)  with  the  estimates  of  the  target  positions  and  intensities. 
The  in-frame  detection  algorithm  realizes  a  generalized  likelihood  ratio  hypothesis  test  with  an  adaptive  threshold 
which  stabilizes  the  false  alarm  rate  and  makes  the  density  of  false  blips  approximately  uniform  over  the  image. 

The  processed  data  from  the  output  of  the  Target  Detection  Block  come  to  the  tracking  system.  The  multitarget 
tracking  scheme  (including  track  management)  is  more  or  less  standard  in  terms  of  the  sequence  of  operations  (see, 
e.g.,[l]).  However,  the  following  important  innovations  are  used:  (1)  global  data  association  (optimal  for  association 
of  all  detections  but  not  locally  optimal  for  a  particular  detection);  (2)  adaptive  selection  of  the  polynomial  power  for 
current  conditions  through  introduction  of  virtual  tracks  in  the  polynomial  filter;  and  (3)  optimal  procedures  for  track 
initiation  and  termination  based  on  changepoint  detection  methods  (see  [2,  3,  5]). 

The  tracker  performs  the  following  operations:  (1)  Identification/association  of  new  detections  (blips)  with  existing 
tracks;  (2)  Initiation  of  new  tracks  based  on  blips  that  are  not  identified  as  belonging  to  existing  tracks;  (3)  Confirmation 
of  newly  initialized  tracks;  (4)  Deletion  of  unconfirmed  tracks;  (5)  Termination  of  tracks;  (6)  Dynamic  correction  of 
the  information  data  base. 


2.3.  Efficiency  Analysis.  We  now  illustrate  the  efficiency  of  the  developed  methodology  by  comparing  with  the 
typical  state-of-practice  differencing  method.  The  differencing  CLUR  method  subtracts  two  consecutive  frames. 

We  simulated  an  image  sequence  with  moderately  intense  clutter  and  sensor  noise  standard  deviation  ctjv  =  3. 
Two  weak  targets  were  inserted  in  the  sequence.  We  first  used  the  Wavelet  spatial-temporal  filter  with  window  of 
m  =  20  frames.  The  results  were  successful  —  the  standard  deviation  of  the  residual  clutter  plus  noise  was  about  3 
and  both  targets  were  tracked,  as  can  be  seen  in  Figure  2(b).  By  contrast,  with  the  differencing  method  we  were  not 
able  to  track  targets,  as  seen  from  Figure  2(c).  Squares  with  no  tracks  attached  represent  instantaneous  detections  part 
of  which  are  false,  while  solid  lines  correspond  to  confirmed  target  tracks. 


(a)  Original  frame 


(b)  Spatial-temporal  Wavelet  CLUR  filter 


(c)  Differencing  CLUR  filter 


Figure  2.  The  results  of  target  tracking  using  spatial-temporal  Wavelet  and  differencing  CLUR  filters. 


3.  CLUTTER  SUPPRESSION  IN  NONSTATIONARY  CONDITIONS 


While  the  clutter  suppression  method  described  in  the  previous  section  is  very  effective  for  geostationary  staring  sen¬ 
sors,  it  does  not  seem  completely  amenable  to  HEO  (e.g.,  SBIRS-HIGH),  LEO  (e.g.,  STSS),  and  other  non-stationary 
moving  sensors.  This  method  should  be  substantially  modified  to  be  efficient  in  non-stationary  conditions.  The  CLUR 
algorithm  should  take  into  account  the  following  specific  features:  (a)  Nonstationarity  related  to  sensor  motion  (not 
only  translational  and  rotational  stabilization  is  needed,  but  also  image  extrapolation,  including  compensation  of  per¬ 
spective  changes);  (b)  A  three-dimensional  cloud  cover  model;  (c)  Complexity  of  images. 

The  developed  CLUR  algorithm  performs  estimation  of  the  projective  transformation  based  on  a  linear  approx¬ 
imation  of  the  image  intensity.  This  results  in  a  system  of  8  nonlinear  equations  solved  by  the  Newton  method. 
For  estimation  of  large  projection  distortions,  a  multi-layer  approach  with  Gaussian  pyramids  is  used.  The  developed 
CLUR  algorithm  covers  very  general  models  with  the  images  that  are  3D-to-2D  projections  of  a  plane  with  an  arbitrary 
orientation  in  3D.  The  details  of  the  algorithm  are  described  below. 


3.1.  Estimation  of  the  Projective  Transformation.  It  is  assumed  that  the  intensity  of  the  image  in  the  vicinity  of 
each  pixel  is  described  by  the  linear  model: 


r/  ,  ,  x  T/  N  dl  ,  dl 

I(x  +  dx,y  +  dy)  =  I{x,y)  +  —dx  +  — 


dy. 


In  order  to  find  a  transformation,  we  minimize  a  “loss”  function  that  is  given  by  the  following  quadratic  norm: 


£2  =  Y  [toihj)  -  I(Px(i,j),Py(i,j))f  =  Y 


i,3 


1,3 


dl  dl 

~  ~  Q^(py(i,j)  -  j) 


(6) 


where  Iq  is  the  reference  frame  and  I  is  the  current  frame. 

Assume  first  that  transformation  leads  to  pixel  shifts  less  than  a  pixel  size,  in  which  case 


Px{i,j) 


cpi  +  ci  j  +  c2 
cpi  +  c7j  +  1  ’ 


Pyihj) 


c%i  +  c^j  +  C5 
cqi  +  c7j  +  1  ' 


(7) 


We  now  have  to  find  the  vector  of  coefficients  c  ==*  (cq  ,  c\ , . . . ,  c7) .  To  this  end,  we  substitute  (7)  in  (6)  and  differentiate 
the  resulting  loss  function  over  coefficients.  The  resulting  derivatives  /o(c),  /i(c), . . . ,  /7(c)  are  equalized  to  zero,  so 
that  we  obtain  the  system  of  eight  nonlinear  equations:  /o(c)  =  0,  /i(c)  =  0, . . . ,  /7(c)  =  0,  where 


/o(c)  =  Y 

i,3 


dl 


dl  . 
dy J 


dl  c0i  +  ci  j  +  c2 
dx  cqI  +  cjj  +  1 


dl  c^i  +  C4J  +  C5  \  dl  i 
dy  c6i  +  Cjj  +  1  )  dx  c6i  +  c7j  +  1  ’ 


For  brevity’s  sake  we  present  only  /o(c).  The  rest  of  the  functions  have  an  essentially  similar  form. 


This  system  is  solved  iteratively  by  the  Newton  method.  Specifically,  let  W(c)  =  \\dfk/dcm\\  be  the  Hessian  of 
the  loss  function  £2  and  let  f(c)  =  (/o(c),  /i(c), . . . ,  /7(c)).  The  ( n  +  l)-st  step  of  the  iterative  method  is  given  by 

^n+l  —  Cn  (Cn)f  (^n) 5  ^  —  0,  1,  ...  . 

Recall  that  we  assumed  that  shifts  are  smaller  that  a  pixel  size,  which  is  not  realistic  for  many  scenarios  of  interest. 
In  order  to  deal  with  shifts  larger  than  a  pixel,  we  use  a  Gaussian  pyramid  with  twice  smaller  images.  The  previous 
procedure  is  used  for  estimation.  The  resulting  matrix  is  scaled  by  two  and  is  then  used  for  image  compensation  at  the 
next  level  of  the  pyramid  with  higher  resolution.  This  process  is  repeated  till  the  original  resolution  is  achieved. 

3.2.  Clutter  Estimation  and  Rejection.  Let  /i, . . . ,  £v  be  the  batch  of  N  sequential  frames  and  let  M^+i  denote 
the  matrix  of  the  projective  transformation  between  the  (i  +  l)-st  and  i- th  frames,  so  that  I*+i  =  M^+iI*  for  i  = 
1, . . . ,  TV  —  1.  Clutter  estimation  is  performed  in  the  coordinate  system  of  the  last  N-th  frame.  To  this  end,  the  matrices 
of  the  projective  transformation  are  multiplied  to  obtain  the  final  transformation  (each-frame-to-the-last- frame): 

M*  =  M iX  Mi+i  x  •  •  •  x  Mat. 

After  all  matrices  are  determined,  clutter  is  estimated  for  each  pixel  in  the  system  of  coordinate  of  the  last  image  by 
the  least  squares  method  (LSM).  For  the  quadratic  model  the  corresponding  system  of  equations  is  linear. 

To  be  specific,  let  the  background  intensity  around  each  pixel  is  approximated  by  the  polynomial  model 

I(dx,  dy )  =  Co  +  C\dx  +  C^dy  +  C%dx2  +  C^dy2  +  C^dxdy. 

where  dx ,  dy  are  coordinates  from  the  center  of  the  output  pixel.  Then  the  LSM  loss  function  is 

N 

C  =  J2(Ii~Co-C1dx- C2dy  -  C3dx2  -  CAdy2  -  C5dxdy ) 2  .  (8) 


More  precisely,  after  all  matrices  are  estimated  the  following  sequence  of  operations  is  performed: 

1.  For  each  point  (u,  v )  we  find  the  closest  point  ( ui ,  Vi)  in  each  frame  in  the  batch  in  the  window  mxm  (usually 
m  =  3),  using  the  estimated  transformation  matrices. 

2.  If  the  number  of  points  is  not  sufficient,  then  the  intensities  are  simply  averaged. 

3.  The  estimate  of  intensity  in  the  current  pixel  is  obtained  from  the  solution  of  linear  equations  (LSM). 

The  third  step  needs  a  little  more  detailed  explanation.  Taking  derivatives  in  (8)  and  equalizing  to  zero,  we  obtain 
the  matrix  equation  for  the  coefficients  Co , . . . ,  C5 .  Note  that  we  are  primarily  interested  in  the  value  of  the  constant 
Co  that  determines  the  intensity  in  the  center  of  the  pixel  in  the  last  frame.  Solving  this  system,  we  find  the  estimate 
Co,  and  then  we  subtract  this  value  from  the  last  frame,  which  yields  the  whitened  frame  (the  output  of  the  clutter 
rejection  filter),  i.e.,  £v  =  In  —  Co- 

3.3.  Extraction  of  Blips  and  Tracking.  For  extraction  of  blips  (detection)  we  use  the  simplest  algorithm  that  per¬ 
forms  thresholding  of  pixel  intensities  in  whitened  images.  The  threshold  is  calculated  based  on  the  following  logic. 
First,  a  histogram  of  intensity  distribution  is  calculated.  Second,  the  histogram  is  integrated  from  small  to  large  in¬ 
tensity  values  until  the  total  number  of  pixels  participating  in  the  integration  does  not  exceed  the  threshold  level. 
For  example,  if  the  frame  size  is  512  x  512  and  the  probability  is  10-4,  then  the  remaining  number  of  pixels  is 
512  x  512  x  10-4  «  26  the  most  intense  pixels.  These  pixels  produce  blips.  This  approach  allows  for  accounting  of 
possible  non-Gaussian  character  of  the  whitened  frame  related  to  residuals  with  clutter  artifacts.  Finally,  the  sub-pixel 
accuracy  improvement  of  the  estimates  of  blip  locations  is  performed.  To  this  end,  in  the  3  x  3  neighborhood  of  each 
blip  the  intensity  is  approximated  with  a  parabola  and  its  minimum  is  found. 

Tracking  is  performed  as  described  in  Section  2.2. 

3.4.  Experimental  Results.  We  performed  a  number  of  simulation  experiments  with  semi- synthetic  data.  In  the 
experiments  targets  are  modeled  as  points  that  are  mixed  with  the  raw  image  prior  to  its  processing.  We  used  IR  image 
sequences  with  real  clutter  (clouds,  ground,  and  a  combination  thereof).  Two  orbits  were  considered:  LEO  and  HEO. 

Figure  3  shows  the  results  of  clutter  rejection  and  target  tracking  for  the  highly  non- stationary  HEO  MeteoSat  data 
with  the  parameters:  Perigee=500km;  Apogee=35800  km;  Input  frame  size  512  x  512.  It  is  seen  that  the  original 
heavy  cloud  clutter  has  been  almost  completely  removed  which  allowed  us  to  detect  and  track  a  dim  point  target  with 


S(C+N)R=0.01.  Green  squares  show  instantaneous  target  detection  in  each  frame  most  of  which  are  false  and  the  red 
line  depicts  the  target  track. 

Figure  4  shows  similar  results  for  the  LEO  scenario  with  the  parameters:  Perigee=200  km,  Apogee=400  km,  Long 
of  Ascending  Node=6.26  rad,  Argument  of  Perigee=0,  Inclination^.  02  rad,  FOV=0.02  rad. 


(a)  Original  frame  with  clouds 


(b)  Whitened  frame  with  target  tracking 


Figure  3.  The  results  of  clutter  suppression  and  target  tracking  for  a  difficult  non- stationary  scenario  (HEO  MeteoSat 
data)  using  the  CLUR  algorithm  with  perspective  prediction  (red  line  -  target  track,  green  boxes  -  blips). 


(a)  Original  frame  with  clouds 


Figure  4.  The  results  of  clutter  suppression  and  target  tracking  for  a  difficult  non- stationary  scenario  (LEO  data)  using 
the  CLUR  algorithm  with  perspective  prediction  (red  line  -  target  track,  green  boxes  -  blips). 


4.  NONLINEAR  FILTERING  BASED  TRACK  BEFORE  DETECT  ALGORITHMS 

In  classical  tracking  methods,  a  thresholding  operation  is  performed  on  the  detection  metric  (e.g.,  based  on  intensity 
or  frame  differences)  formed  in  individual  scans/frames  of  data  to  generate  target  detections  that  are  subsequently  pro¬ 
cessed  by  data  association  and  state  estimation  algorithms.  Thresholding,  however,  causes  partial  loss  of  information, 
and  target  tracking  based  on  this  approach  may  be  completely  infeasible  in  situations  with  severe  clutter  and/or  noise 
(if  the  detection  threshold  is  set  low  enough  to  detect  very  dim  targets,  the  number  of  clutter  detections  may  result  in 
either  poor  track  estimation  performance,  or  may  prohibit  real-time  operation  as  the  number  of  data  association  hy¬ 
potheses  explodes).  In  track-before-detect,  information  from  individual  frames  is  integrated  over  time  until  the  number 
and  locations  of  targets  can  be  accurately  estimated  [2,  6,  7].  In  our  approach  to  track-before-detect,  clutter-suppressed 
frames  are  used  to  update  a  non-parametric  statistical  background  scene  model  that  is  built  up  on-the-fly  to  account  for 
intensity  variations  due  to  effects  such  as  noise  and  clutter  leakage.  Each  clutter-suppressed  frame  is  evaluated  in  the 


statistical  model  to  estimate  the  probability  that  the  intensity  of  each  optical  sample  has  been  modified  by  the  presence 
of  a  moving  target.  The  moving  target  probabilities  are  then  integrated  over  time,  based  on  probabilistic  constraints  on 
possible  target  motions,  to  develop  the  a  posteriori  probability  distribution  function  over  the  space  of  possible  target 
locations,  based  on  the  measurement  history  and  any  available  a  priori  information  (note  that  this  a  priori  information 
is  not  a  requirement).  A  track  can  be  initiated  by  an  external  cue,  such  as  a  track  formed  from  another  sensor  processor, 
or  initiated  automatically  by  our  multiple  target  hypothesis  management  algorithm.  In  the  latter  case,  detections/tracks 
are  declared  and  reported  when  likelihoods  of  hypothesized  target  trajectories  become  sufficiently  large.  Subsequently, 
the  track  state  estimates  are  updated  with  information  from  each  new  data  scan  as  it  becomes  available. 

Estimation  of  the  posterior  probability  over  the  space  of  possible  target  trajectories  is  complicated  by  the  nonlin¬ 
earities  in  target  dynamics  and  optical  measurement  (which  involves  a  highly  nonlinear  projective  function).  To  meet 
this  challenge,  we  employ  particle  filtering  methods — in  particular,  the  sampling-importance-resampling  (SIR)  filter, 
for  modeling  and  updating  the  posterior  distribution  [8].  Particle  filters  are  well-suited  for  this  problem  since  they  are 
essentially  non-parametric  density  estimators,  and  there  is  no  constraint  on  the  underlying  structure  of  the  probability 
density  function  (PDF).  In  contrast,  state  estimation  problems  in  which  nonlinear  dynamic  models,  nonlinear  mea¬ 
surement  models,  or  non-Gaussian  noise  processes  are  present  pose  fundamental  difficulties  when  applying  standard 
Kalman  filtering  techniques  based  on  linearization.  Particle  filtering  algorithms  do  not  require  linearization  of  dy¬ 
namic  equations  or  measurement  equations.  Updating  the  target  state  PDF  simply  requires  values  of  the  measurement 
likelihood  function  evaluated  at  discrete  points  in  the  state  space.  In  our  algorithm,  the  likelihoods  are  obtained  from 
evaluating  each  clutter-suppressed  optical  sample  in  the  statistical  background  model,  as  previously  described.  In  the 
case  of  tracking  in  3D  coordinates,  we  apply  the  sensor  measurement  model  to  project  the  3D  target  position  hypoth¬ 
esis  for  each  particle  into  the  image  plane,  prior  to  likelihood  updating.  For  propagation  through  the  state  dynamic 
model,  each  particle  is  simply  substituted  into  the  dynamic  equation  and,  with  the  simulation  of  any  process  noise 
variables,  a  new  particle  results. 

As  shown  in  the  example  in  Figure  5,  after  multiple  frames  the  particle  distribution  converges  such  that  clusters 
of  particles  are  located  in  the  most  likely  regions  of  the  state  space,  with  a  small  fraction  of  the  particles  distributed 
more  diffusely.  As  a  result,  we  have  found  that  our  particle  filter-based  track-before-detect  algorithm  is  more  than  two 
orders  of  magnitude  more  efficient  (in  terms  of  computation  and  storage  requirements)  compared  with  the  classical 
3D  velocity-matched  filter  [7],  for  the  same  level  (near-optimal  effective  SNR  increase)  of  performance.  This  is  true 
even  though  the  particle  filter  was  configured  to  account  for  possibly  large  target  maneuver  accelerations.  Notably, 
our  current  C++  implementation  of  particle  filtering-based  track-before-detect  is  able  to  process  MSX  SBV  data  in 
near-real-time.  The  explanation  for  the  efficiency  gain  in  the  recursive  particle  filter  algorithm  is  that  after  processing 
each  frame  of  data,  during  the  particle  re-sampling  step,  the  particles  are  automatically  re-allocated  to  perform  more 
finely-spaced  searches  for  target  locations  and  velocities  in  regions  of  the  state  space  that  are  found  to  be  most  likely 
(based  on  processing  multiple  past  frames  of  data),  whereas  in  the  3D  velocity-matched  filter,  hypotheses  are  allocated 
uniformly  while  processing  the  entire  batch  of  frames. 

In  Figure  5,  we  provide  an  example  of  track-before-detect  using  MSX  SBV  data.  This  example  is  interesting 
because  it  contains  two  resident  space  objects  (RSOs),  one  which  is  blurred  by  rapid  motion  of  approximately  30 
pixels/frame,  and  the  other,  which  is  much  dimmer,  moves  at  only  approximately  1/3  pixel/frame.  This  requires 
a  very  diffuse  distribution  of  particle  velocities,  and  fairly  large  process  noise  power  spectral  density  in  the  target 
dynamics  model.  Furthermore,  it  is  important  to  be  able  to  track  both  the  dim  (mean  intensity  of  3)  and  bright  (mean 
intensity  of  19)  targets  simultaneously,  which  is  possible  using  our  statistical  background  modeling  algorithm,  which 
accurately  recognizes  both  low-  and  high-intensity  moving  targets  as  anomalies,  compared  with  the  background  scene. 
In  Fig.  5(a),  the  fifth  frame  in  the  sequence  is  shown  (the  square  root  of  the  intensities  are  plotted  to  provide  improved 
contrast).  In  Fig.  5(b),  the  initial  uniform  distribution  of  particle  locations  is  shown  (in  this  example,  a  2D  linear 
velocity  white  noise  acceleration  target  dynamics  model  was  used),  with  lighter  intensities  denoting  larger  numbers  of 
particles  within  an  optical  resolution  cell.  In  Fig.  5(c),  the  target  probabilities  produced  by  evaluating  the  fifth  frame  in 
the  statistical  background  model  are  shown,  with  lighter  intensities  denoting  higher  probabilities.  Finally,  in  Fig.  5(d), 
the  particle  distribution  after  processing  the  fifth  frame  is  shown,  and  cues  are  provided  to  the  locations  of  both  targets, 
which  have  been  automatically  detected  (after  3  and  5  frames  for  the  brighter  and  dimmer  targets,  respectively). 

In  Figure  6,  we  provide  an  example  of  track-before-detect  in  3D  coordinates.  The  frame  sequence  was  generated 
synthetically  using  measured  EUMETSAT  clutter,  but  with  simulated  target  trajectories  and  intensities  (SCR  =  1/1000) 
and  additive  white  Gaussian  sensor  noise  (SNR  =  1/6).  As  shown  in  this  example,  tracking  in  3D  physical  world 
coordinates  using  a  single  optical  sensor  is  possible  in  some  cases  (e.g.,  non-nadir  viewing  of  a  ballistic  trajectory  with 
a  known  launch  point),  and  when  ambiguity  does  exist,  it  is  well-represented  in  the  particle  filter.  In  this  example, 


(a)  Fifth  frame  (square  root  contrast  enhanced) 


(c)  Moving  target  probabilities  (for  fifth  frame) 


Figure  5.  Track-before-detect  results  in  MSX  SBV  for  two  RSOs  (one  extremely  dim)  with  velocities  and  intensities 
differing  by  factors  of  100  and  6,  respectively. 


we  assumed  that  the  launch  location  and  time  were  known  (e.g.,  from  bright  plume  detection),  and  that  an  initial 
estimate  of  the  missile  launch  heading  and  elevation  angles  was  available.  The  angle  estimates  were  assumed  to 
be  very  inaccurate,  such  that  the  initial  uncertainty  on  the  impact  point  location  was  very  large:  1,500  km  3-sigma 
uncertainty,  for  a  trajectory  length  of  1,500  km.  In  Fig.  6(b),  the  particle  distribution  is  shown  after  processing  50 
frames  at  a  rate  of  1  Hz,  and  the  presence  of  a  single  target  is  clearly  seen  in  the  peaked  distribution.  The  target 
has  been  automatically  confirmed  by  the  algorithm  before  this  time,  although  visual  cueing  was  not  used  in  this 
experiment.  In  general,  tracking  in  3D  provides  improved  estimation  accuracy,  and  enables  improved  multi-sensor 
fusion  and  trajectory  prediction,  and  as  shown  here,  when  combined  with  sensor  motion  compensation  and  spatial- 
temporal  clutter  suppression,  this  processing  can  be  applied  even  for  dim  targets  against  terrestrial  and  other  highly 
cluttered  backgrounds. 


5.  CONCLUSIONS 

In  this  paper,  we  have  shown  the  need  for  and  effectiveness  of  data-driven  registration  processing,  spatial-temporal 
clutter  estimation  and  rejection,  and  nonlinear  stochastic  track-before-detect  processing  for  low  observable  targets 
against  cluttered  terrestrial,  earth-limb,  and  celestial  backgrounds.  Notably,  we  have  demonstrated  that  for  some  real- 
world  applications  it  is  currently  feasible  to  perform  this  processing  in  near-real-time  using  a  desktop  PC.  Furthermore, 


(a)  Simulated  frame  with  EUMETSAT  background  (b)  Projection  of  3D  particles  into  image  plane 

clutter 


Figure  6.  Projection  of  3D  particles  (after  50  sec.)  into  the  image  plane  for  a  simulated  target  against  EUMETSAT 
background  clutter  (SCR  =  1/1000)  with  simulated  sensor  noise  (SNR  =  1/6). 


we  expect  that  the  developed  algorithms  have  utility  in  a  wide  range  of  current  and  future  systems  being  developed  for 
space-based  EO  and  IR  surveillance. 
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